Time Series Analysis
Functions and Examples
By Majid M Al-Sadoon (majid.sadoon@aramco.com.sa) 11/3/2003
Based on "Time Series: Theory and Methods" 2nd Ed. by Brockwell and Davis

Sample Autocovariance p.29

\(g[V_, h_] := 
      Module[{m, n}, 
		n = \(Dimensions[V]\)[\([1]\)]; 
		m = 
          Sum[\(V[\([t]\)]\)[\([2]\)], {t, 1, n}]/
            n; 
		Sum[\((\(V[\([t + Abs[h]]\)]\)[\([2]\)] - 
                  m)\) \((\(V[\([t]\)]\)[\([2]\)] - m)\), {t, 1, n - Abs[h]}]/
          n]\)

Yule-Walker Experiment - calls g p.239

<< Statistics`ContinuousDistributions`

ndist = NormalDistribution[0, 2]

w = Table[{i, Random[ndist]}, {i, 1, 200}];

z[i_] := w[[i]][[2]]
x[i_] := -0.8 x[i - 1] + z[i]
\(	\(x[1] = 1;\)\)

X = Table[{i, x[i]}, {i, 1, 200}];

ListPlot[X, PlotJoined -> True]

g1 = Table[{h, g[X, h]/g[X, 0]}, {h, 0, 40}];
ListPlot[g1, PlotJoined -> True]

Γ1 = 
        Table[g[X, i + j - 2], {i, 1, 1}, {j, 1, 1}];
γ1 = Table[g[X, i], {i, 1, 1}];

Inverse[Γ1] . γ1

Innovations Algorithm - calls g p.245

Innovations[V_, m_] := 
      Module[{v, x, θ, H, A}, 
x[0] = 
          0; 
v[0] = 
          g[V, 0]; 

  
θ[n_, 
            l_] := θ[n, l] = 
            If[l == 0, 
              1, (g[V, l] - 
                    Sum[θ[n - l, n - l - j]*θ[n, n - j]*
                        v[j], {j, 0, n - l - 1}])/
                v[n - l]]; 

v[
            n_] := v[n] = 
            v[0] - Sum[(θ[n, n - j]^2)*v[j], {j, 0, 
                  n - 1}]; 

x[
            n_] := x[n] = 
            Sum[θ[n, 
                  j] (V[[n + 1 - j]][[2]] - x[n - j]), {j, 1, 
                n}]; 

H[n_] := 
          Table[θ[n, i], {i, 1, 
              n}]; 

A[n_] := 
          Table[Sum[θ[n, i - r] θ[n, j - r], {r, 1, 
                Min[i, j]}], {i, 1, n}, {j, 1, 
              n}]; 

Table[{x[n], v[n], 
            H[n], A[n]/Dimensions[V][[1]]}, {n, 0, 
            m}]
]

Needs["Statistics`ContinuousDistributions`"]

ndist := Random[NormalDistribution[0, Sqrt[2.25]]]

hh = {-1.40, 0.60};
Solve[z^2 + hh . {z, 1} == 0, z]
W = Table[ndist, {i, 1, 500}];
ClearAll[y];
y[0] = 0;
y[t_] := 
        W[[t + 2]] + 
          Evaluate[{W[[t + 1]], W[[t]]} . 
              hh];
V = Table[{t, y[t]}, {t, 0, 498}];
ListPlot[V, PlotJoined -> True, 
      AxesOrigin -> {1, 0}]
ListPlot[Table[{h, g[V, h]/g[V, 0]}, {h, 0, 10}], 
      PlotJoined -> True, PlotRange -> {-1, 1}]

Innovations[V, 5] // TableForm

Preliminary Estimates for an ARAM (p,q) - calls Innovations p.250

PreliminaryARMA[V_, p_, q_, n_] := 
      Module[{H, hm, h, f, eq}, 
If[p + q > n, 
          Print["Error! p + q > ", 
            n], 
H = Innovations[V, 
                  n][[n + 1]][[3]]; 
hm[0] = 
            1; 
hm[i_] := 
            H[[i]]; 
h[0] = f[0] = 
              1; 
eq = 
            Table[hm[j] == 
                If[j ≤ q, h[j], 0] + 
                  Sum[If[i ≤ p, f[i], 0]*hm[j - i], {i, 1, 
                      Min[j, p]}], {j, 1, 
                p + q}]; 
{Flatten[{Table["f"[j], {j, 
                    1, p}], Table["h"[j], {j, 1, q}]}], 
            Flatten[{Table[f[j], {j, 1, p}], 
                  Table[h[j], {j, 1, q}]}] /. Solve[eq, 
                  Flatten[{Table[f[j], {j, 1, p}], 
                      Table[h[j], {j, 1, 
                          q}]}]][[1]]}
]]

Needs["Statistics`ContinuousDistributions`"]

ndist := Random[NormalDistribution[0, Sqrt[2.25]]]

W = Table[ndist, {i, 1, 500}];
ClearAll[y];
y[0] = -1.1;
y[t_] := 
        W[[t + 2]] - 
          1.4  W[[t + 1]] +  .6  W[[t]];
V = Table[{t, y[t]}, {t, 0, 98}];

PreliminaryARMA[V, 0, 2, 30] // TableForm

Durbin-Levinson Algorithm - calls g p.242

DL[V_, m_] := 
      Module[{x, ϕ, 
          v, Φ, Γ}, 
x[0] = 
          0; 
ϕ[1, 1] = 
          g[V, 1]/g[V, 0]; 
v[1] = 
          g[V, 0] (1 - ϕ[1, 
                    1]^2); 

ϕ[n_, 
            j_] := ϕ[n, j] = 
            If[n == 
                j, (g[V, n] - 
                    Sum[ϕ[n - 1, i] g[V, n - i], {i, 1, n - 1}])/
                v[n - 1], ϕ[n - 1, 
                  j] - ϕ[n, n] ϕ[n - 1, 
                    n - j]]; 

v[
            n_] := v[n] = 
            v[n - 1] (1 - ϕ[n, 
                      n]^2); 

x[
            n_] := x[n] = 
            Sum[ϕ[n, j] V[[n - j + 1]][[2]], {j, n, 
                1, -1}]; 

Φ[n_] := 
          Table[ϕ[n, j], {j, 1, 
              n}]; 

Γ[n_] := 
          Table[g[V, i - j], {i, 1, n}, {j, 1, 
              n}]; 

Table[{x[n], 
            v[n], Φ[n], 
            v[n] Inverse[Γ[n]]/Dimensions[V][[1]]}, {n, 
            1, m}]]

Needs["Statistics`ContinuousDistributions`"]

ndist := 
        Random[NormalDistribution[0,  .2]];
W = Table[ndist, {i, 1, 100}];
Clear[y];
y[0] = 0;
y[1] = 0;
y[
          t_] := y[t] =  .1  y[t - 1] -  .1  y[t - 2] + 
            W[[t]];
V = Table[{t, y[t]}, {t, 0, 99}];
ListPlot[V, PlotJoined -> True, 
      AxesOrigin -> {1, 0}]
ListPlot[Table[{h, g[V, h]/g[V, 0]}, {h, 0, 15}], 
      PlotJoined -> True, PlotRange -> {-1, 1}]

DL[V, 4] // TableForm

Sample ACV and PACV - calls DL p.243

ACV[V_, m_] := g[V, m]/g[V, 0]

PACV[V_, m_] := 
      If[m == 0, 
        1, DL[V, m][[m]][[3]][[m]]]

Calculation of ψ(z) = θ(z) / ϕ(z) p.91

Ψ[ϕ_, θ_, n_] := 
      Module[{ψ, f, h}, 
f[0] = 
          1; 
f[i_] := f[i] = 
            If[i ≤ Dimensions[ϕ][[1]], ϕ[[i]], 0]; 
h[0] = 
          1; 
h[i_] := h[i] = 
            If[i ≤ Dimensions[θ][[1]], θ[[i]], 0]; 
ψ[j_] := ψ[j] = 
            If[0 ≤ j < 
                Max[Dimensions[ϕ][[1]], 
                  1 + Dimensions[θ][[1]]], 
              h[j] + Sum[f[k] ψ[j - k], {k, 1, j}], 
              Sum[f[k] ψ[j - k], {k, 
                  1, Dimensions[ϕ][[1]]}]]; 
Table[ψ[i], {i, 0, n}]]

Ψ[{ .25}, {}, 10] // N
Ψ[{1, - .25}, {1}, 10] // N

ψ[
          z_] := (1 + z)/(1 - 
              z +  .25  z^2);
InverseZTransform[ψ[1/z], z, 4]

Autocovariance of an ARAM (p,q) - calls Ψ p.93

γ[ϕ_, θ_, σ_, n_] := 
      Module[{eq, h, K, G, ψ, gam}, 
h[0] = 
          1; 
h[i_] := h[i] = 
            If[i ≤ Dimensions[θ][[1]], θ[[i]], 
              0]; 
K = Ψ[ϕ, θ, 
            Max[Dimensions[ϕ][[1]], 
              1 + Dimensions[θ][[1]]]]; 
ψ[j_] := K[[j + 1]]; 
eq = 
          Table[{Γ[k] - 
                  Sum[ϕ[[j]] Γ[Abs[k - j]], {j, 
                      1, Dimensions[ϕ][[1]]}] == (σ^2) Sum[h[j] ψ[j - k], {j, 
                      k, Dimensions[θ][[1]]}]}, {k, 
              0, Dimensions[ϕ][[1]]}]; 
G = (Table[Γ[k], {k, 0, Dimensions[ϕ][[1]]}] /. 
                Solve[Flatten[eq], 
                  Table[Γ[k], {k, 
                      0, Dimensions[ϕ][[1]]}]])[[1]]; 
        For[k = 0, k <= Dimensions[ϕ][[1]], 
          gam[k] = G[[k + 1]]; k++]; 
gam[
            k_] := gam[k] = 
            Sum[ϕ[[j]] gam[k - j], {j, 
                  1, Dimensions[ϕ][[1]]}] + 
If[
                0 <= k < 
                  Max[Dimensions[ϕ][[1]], Dimensions[θ][[1]] + 1], σ^2*
                  Sum[h[j] ψ[j - k], {j, 
                      k, Dimensions[θ][[1]]}], 
                0]; 
Table[gam[j], {j, 0, n}]]

γ[{}, {θ}, σ, 3]

γ[{ϕ}, {θ}, σ, 3] // Simplify

γ[{ϕ}, {}, σ, 3]